function [foE, NmE] = foE_NmE(lat, Az, chi_eff, mth)
%
% E-Layer Critical Frequency
%
%DESCRIPTION:
%This function computes the E-Layer Critical Frequency (foE, [MHz]) at a
%given location, in addition to the effective Solar Zenith Angle (chi_eff),
%a season dependent parameter has to be computed.
%
%PROTOTYPE:
% [foE, NmE] = foE_NmE(lat, Az, chi_eff, mth)
%
%--------------------------------------------------------------------------
% INPUTS:
%   lat        [1x1]       Latitude                  [deg]
%   Az         [1x1]       Eff. Ionisation Level     [-]
%   chi_eff    [1x1]       Eff. Solar Zen. Ang.      [deg]
%   mth        [1x1]       Month                     [month]
%--------------------------------------------------------------------------
% OUTPUTS:
%   foE        [1x1]       E-Layer Crit. Freq.       [MHz]
%   NmE        [1x1]       E-Layer Max. Density      [1e+11 m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Constants
DR = pi/180;                %Conversion Factor: deg->rad

%Definition of the Seas Parameter (as func. of month of the year)
if     (mth == 1) || (mth == 2) || (mth == 11) || (mth == 12)
    seas = -1;
elseif (mth == 3) || (mth == 4) || (mth == 9)  || (mth == 10)
    seas = 0;
elseif (mth == 5) || (mth == 6) || (mth == 7)  || (mth == 8)
    seas = 1;
end

%Introduction of Latitude Dependence
ee = exp(0.3*lat);
seasp = seas * (ee-1)/(ee+1);

%E-Layer Critical Frequency [MHz]
foE = sqrt((1.112 - 0.019*seasp)^2*sqrt(Az)*(cos(chi_eff*DR))^(0.6) + 0.49);

%E-Layer Maximum Density [10^11 m-3]
NmE = 0.124 * foE^2;

end